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Abstract 

The  search  for  a  universal  solution  of  the  equations  of  motion  for  a 
satellite  orbiting  an  oblate  planet  is  a  subject  that  has  merited  great 
interest  because  of  its  theoretical  and  practical  implications.  Here,  a 
complete  first-order  perturbation  solution,  including  the  effects  of  the 
J2  terms  in  the  planet's  potential,  is  given  in  terms  of  standard  orbital 
parameters.  The  simple  formulas  provide  a  fast  method  for  predict¬ 
ing  satellite  orbits  that  is  more  accurate  than  the  two-body  formulas. 

These  predictions  are  shown  to  agree  well  with  those  of  a  completely 
numerical  code  and  with  actual  satellite  data.  Aiso.  in  an  appendix,  it 
is  rigorously  proven  that  a  satellite  having  negative  mechanical  energy 
remains  for  all  time  within  a  spherical  annulus  with  radii  approximately 
equal  to  the  perigee  and  apogee  of  its  initial  osculating  ellipse. 
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1  Introduction 


A  characteristic  feature  of  practical  orbit  prediction  is  that  the  engineer  may  deal  with 
numerous  satellites  in  a  great  variety  of  orbits.  Under  these  circumstances  analytical  relations 
which  can  quickly  approximate  an  orbit  may  be  far  superior  to  large  numerical  programs. 
While  many  analytical  models  have  been  developed  for  the  artificial  satellite  age.  most  are 
not  used  in  practical  orbit  prediction  because  they  violate  one  or  more  of  the  following 
principles: 

•  The  method  should  provide  a  solution  that  is  significantly  more  accurate  than  the 
two-body  solution. 

•  The  real  physical  effects  of  the  orbit  should  be  easily  distinguishable  in  the  solution. 

•  The  solution  should  be  universal:  it  should  be  valid  for  all  orbital  parameters. 

The  problem  of  predicting  the  motion  of  a  satellite  perturbed  only  by  the  oblateness  of  the 
planet  has  received  considerable  attention  following  the  first  launchings  of  artificial  satellites 
about  the  Earth.  Some  of  the  studies  of  this  problem  by  means  of  general  perturbation 
theories  are  listed  at  the  end  of  this  paper.  Techniques  have  involved  expansions  in  powers 
of  y/T2.  averaging  processes,  the  use  of  spheroidal  coordinates,  and  the  edifice  of  Hamiltonian 
mechanics.  It  is  not  the  intention  of  this  present  paper  to  compare  the  various  methodologies 
used.  Suffice  it  to  say  that  many  researchers  believe  a  solution  which  embodies  all  of  the 
above  principles  was  not  achieved  (e.g..  see  Taff). 

The  basic  procedure  used  in  this  paper  to  solve  the  differential  equations  of  motion  is 
the  perturbation  technique  known  as  the  Method  of  Strained  Coordinates.  This  technique 
was  first  applied  to  the  title  problem  by  Brenner.  Latta.  and  Weisfield.  Using  a  rruan  orbital 
plane  to  specify  an  arbitrary  orbit,  they  were  only  able  to  obtain  a  partial  solution  (e.g..  the 
eccentricity  was  assumed  small  and  initial  conditions  were  not  considered). 


Here  we  use  coordinates  in  the  true  orbital  plane  to  cast  the  differential  equations  into  a 
simplified  form,  as  was  originally  done  by  Struble. 

2  Orbital  Kinematics 

Figure  1  shows  the  usual  reference  system  of  spherical  coordinates  (r, o.i?).  The  radial 
distance  r  is  measured  from  the  center  of  the  planet  0  to  the  satellite  5.  The  line  is  in 
a  direction  fixed  with  respect  to  an  inertial  coordinate  system.  The  right  ascension  a  is  the 
angle  measured  in  the  planet's  equatorial  plane  eastward  from  the  line  0~j.  The  declination 
or  latitude  j3  is  the  angle  measured  northward  from  the  equator.  The  position  vector  r  of 
the  satellite  in  the  spherical  coordinate  system  is 

r  =  r(cosa  cos  J)bj  -f  r(sin  a  cos  i)b2  +  r(sin  i)b3  ( 1 ) 

where  (bj.b2.b3)  are  orthonormal  base  vectors  fixed  in  the  directions  shown. 

We  can  also  locate  the  satellite  by  its  polar  coordinates  (r.  6)  within  a  (possibly  rotating) 
orbital  plane  that  instantaneously  contains  its  position  and  velocity  vectors.  Here  0  is  the 
argument  of  latitude,  i.e.,  the  angle  measured  in  the  orbital  plane  from  the  ascending  node  to 
the  satellite.  The  orbital  plane  is  inclined  at  an  angle  i  to  the  equatorial  plane  and  intersects 
the  equatorial  plane  in  the  line  of  nodes,  making  an  angle  Q  with  the  O'  line. 

We  introduce  another  orthonormal  set  of  base  vectors  (Bj.B2.B3)  which  move  with  the 
satellite  so  that  Bj  is  in  the  direction  of  the  position  vector  r.  B2  is  also  in  the  orbital  plane, 
and  B3  =  Bj  x  B2.  The  basis  (bj .  b2.  b3)  may  be  transformed  into  the  basis  (Bj .  B2.  B3)  by 
a  succession  of  three  rotations.  First  the  basis  (bj,b2,b3)  is  rotated  about  the  b3  direction 
by  the  angle  SI,  next  the  basis  is  rotated  about  the  new  1-direction  by  the  angle  i.  and 
finally  the  basis  is  again  rotated  about  the  new  3-direction  by  the  angle  0.  The  two  sets  of 
base  vectors  are  related  by  the  product  of  the  rotation  matrices  representing  each  successive 
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rotation  (as  explained  in  the  book  by  Danielson): 
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The  position  vector  r  has  only  one  component  in  the  rotating  basis: 

r  =  rBj  (3) 

Using  the  first  of  equations  (2).  we  obtain  the  components  of  r  in  the  fixed  basis: 

r  =  r(cos  #cos  Cl  —  sin  9  cos  i  sin  Cl ) bx 

-fi  r(cos  9  sin  Q  4-  sin  9  cos  ?  cos  fljbj  +  r(sin  9  sin  ?  )b3  (4 ) 

Equating  the  components  of  equations  (1)  and  (4).  we  can  obtain  the  following  relations 
among  the  angles  (q.3)  of  the  spherical  coordinate  system  and  the  agronomical  angles 

{i.n.ey. 


sin  3  =  sin  #sin  t 


(o) 


cos  3  =  cos  9  sec(a  -  Cl) 


The  velocity  dr/dt  of  the  satellite  is  obtained  by  differentiating  (3)  with  respect  to  the 
time  t : 

dr  dr  _  dB, 

(6) 


dr  dr  dB] 

W=d?B'  +  r7T 


Since  the  orbital  plane  must  contain  the  velocity  vector,  we  have  to  enforce 

dBj 


dt 


b3  =  o 


(<) 


Substitution  of  equations  (2)  into  equation  (7)  leads  to  a  relationship  which  uncouples  the 
equations  for  Cl(9)  and  ?(9): 

dCl  tan  9  di 

d$  sin  i  d0  ' 
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The  velocity  (6)  can  then  be  written 


dr  dr  d6(  di  \ 

—  =  — Bi  4-  r—  1  4-  tan  6  cot  ?  —  Bj  (9) 

dt  dt  dt  \  ddj 

In  the  following  part  of  this  paper,  we  will  obtain  expressions  for  r(8 j.  ?(#).  fl(8).  and 
dt/d0(9).  The  position  and  velocity  vectors  of  the  satellite  then  may  be  calculated  from 
the  formulas  in  this  section.  The  classical  orbital  elements  p.t.  and  u.'  are  the  semilatus 
rectum,  eccentricity,  and  argument  of  perigee  of  the  instantaneous  (osculating)  conic  section 
determined  by  the  position  and  velocity  vectors.  If  needed.  p($).  e(0).  and  can  be 

obtained  from  our  solution  r(6)  and  dtjd9{6\. 


(  cos(6  - 


(  sin(#  -  _■ 


.A  -  L  it 


3  Equations  of  Motion 


The  expressions  for  the  kinetic  and  potential  energies  per  unit  mass  of  a  satellite  orbiting 


around  an  oblate  planet  are  respectively: 


T  1  (drW 

T=2  U  +r 


2  f -r-  V  +  r2 cos2  J  (~\2 


GM  J2R‘'  .  2 

l  = — —  [l  +  ^T  (2  ~  3sm  3) 


where  G  is  the  gravitational  constant.  M  is  the  mass  of  the  planet,  R  is  the  equatorial  radius 
of  the  planet,  and  J2  is  the  constant  coefficient  of  the  spherical  harmonic  of  degree  2  and 
order  0  in  the  planet's  gravitational  field.  Substitution  of  these  equations  into  Lagrange's 


equations 


d  d(T-V)  d_ 

dt  a(£) 


(T  -  n  =  o 


q  =  r.  a.  or  3 


results  in  the  following  equations  of  motion; 


(d3\2  _ 

U )  r 


cos2  3 


j-  (r2  cos2  3~)  =  0 
dt  l  dt 


)  -f  r2  sin  5  cos  3  (  -77  =  ~ 


dt  \  -  /  93 

Initial  conditions  are  established  by  requiring  that  at  the  initial  time  t0  the  orbital  pa¬ 
rameters  of  the  usual  two- body  orbit,  the  conic  section  determined  by  the  initial  position  and 
velocity  vectors,  are  known.  The  actual  orbit  is  then  tangent  to  this  initial  instantaneous 
conic  section  at  to  (see  Figure  1).  Equating  the  initial  position  and  velocity  vectors  given  by 
equations  (3)  and  (9)  to  the  two-body  expressions,  we  obtain 


r(t  0)  = 


1  +  e0cos(&0  - 


dr  x  f  0^0  sinffly  -  w0) 

-  - 

dt  p0 


d9,  , 


ro  +  tan  cot  t0^(^o)| 
i($o)  =  >0 


n(o0)  =  o0  as) 

Here  ho  =  \jGMpo  is  the  initial  value  of  the  satellite’s  specific  angular  momentum  about  the 
center  of  the  planet,  and  the  subscript  0  on  a  symbol  denotes  that  the  parameter  is  evaluated 
at  the  initial  time  t0- 

We  immediately  have  two  integrals  of  the  equations  of  motion: 


T  +  V  =  constant 


2  2  -  dci 

r  cos  3~—  —  constant 
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Equation  (19)  simply  states  that  the  mechanical  energy  of  the  satellite  remains  constant. 
Now,  from  equations  (1)  and  .16) 


* 


2  2  Jo  dr 

r  cos  +—  =  r  x  —  •  b3  =  h0  cos  i0 
at  at 


(21) 


Equation  (21)  simply  states  that  the  component  along  the  polar  axis  of  the  specific  angulai 
momentum  of  the  satellite  remains  constant.  Inserting  equations  (3)  and  (9)  into  equation 
(21 ),  we  obtain 


d t_ 
dd 


r 2  cos i 
ho  cos !0 


d? 

1  +  tan  0  cot  i  — 
du 


( 


•70 


! 


This  allows  the  independent  variable  to  be  changed  from  1  to  0. 

Letting  u  —  Po/r-  and  using  equations  (5).  (21 .).  and  (22).  we  can  rewrite  the  remain  ina 
equations  of  motion  (12)— (13): 


di  —  2Jv  sin  0  cos  0  sin  i  cos2  t 
d&  si — f  2Ju  sin2  0  cos3  i 

COS  J 


d2u  cos2 1  Jcos2i 

-  “f*  1/  —  -  -j-  — — — * - 

dB  2  C2  C2 


t.2(  1  —  3  sin*  6  sin2  i)  +  2 u~  sin  0  cos  8, 1  —  3  cos2 1  j 

dB 


sin2  0cos2  i  -  2  (  sin2  0  cos2; 
du-  \  du  I 


4  J2v  sin3  9  cosf'  i 


du  7 

it  —  cos 0(2  +  sin  i) 
dd 


d  (  du 

1 )  +  17,  1  v~  '  s,n  "  cos*  ’ 


de  \  de 


The  terms  in  (24)  with  d2u/d9 2  can  be  combined,  yielding  the  equivalent  equation 

d2u  f  cos2  i  J  cos2  i .  5  .  ,  „  ,  . 

•7  +  ti=  — r - 1- - 5 — u2(l  +  sin*  0{  t  cos*  i  -  3)) 

du-  1  cz  c 2 
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du  (  du.  \  ^ 

+2u— —  sin  B  cos  0(1  —  3  cos2 1)  —  2  (  —  )  sin2  0cos2 1) 
dO  \ aB J 


+ 


4  J2u  sin3  9  cos6  ? 


u2  sin  0  cos*  i  —  u  cos  0(2  4-  sin2  ?) 

UP 


du 

dd 


sin  0  cos2  t 


4(1  + 


4  Ju  sin2  0  cos4 ;  4  J2u2  sin4  0  coss 


•  ? 


i  25) 


Here  we  have  introduced  the  shorthand  notation  c  =  cosio,  $  =  sin  i0.  J  —  3J2/?2/2/)2. 


4  Perturbation  Procedure 


The  differential  equations  (23 )— •  24)  are  coupled  by  the  nonlinear  terms  and  apparently 
cannot  be  solved  analytically.  If  we  expand  the  right  sides  of  (23)  and  (25)  in  a  Taylor  series 
expansion  in  powers  of  J  and  retain  only  terms  up  to  order  ./*.  the  equations  simplify  to 


di  —2Ju  sin  6  cos  6  sin  ;  cos3  4  J2u2  sin  ?  cos' 


-  sin3  9  cos  9  4-  0(  JA 


+  V  = 


cos2;  J  cos2 ;  f  —  4u  sin2  9  cos4 


4  u2[l  4-  sin2  9(  7  cos2  ;  —  3); 


du  .  .  „  -j  ( du\  ,  ,  , 

4-  2u  —  sin  9  cos  6(1  —  3  cos2 /  )  —  2  —  )  sin'  9  cos2  ; 
du  \  du  ! 


4  J2u  sin2  9  cos6 


u2\~  1  +  3  sin2  0(  1  -  2  cos2  i )]  + 


3u  sin2  9  cos4 


du  .  „  „ ,  t  .  ,  /  du  \  *  .  t  ^ 

+  u  —  sin  6  cos  0.  i  cos* ;  -  5j  4  I  —  1  sin*  6  cos" ;  >  4-  0(d  ) 

Here  the  term  in  the  0  svmbols  indicates  that,  for  all  sufficiently  small  J.  the  err^'  is 


is  less 


than  a  constant  times  Jz .  The  equations  (26)-(27)  are  identical  to  those  used  as  the  starting 
point  in  the  analysis  of  Eckstein,  et  al. 

It  is  reasonable  to  expect  that  the  solution  for  u  will  be  arbitrarily  close  to  tin1  tvo  body 
solution.  1  4-  e0  cos ( ^  —  -Vo),  when  J  is  close  to  zero.  This  assumption  is  consistent  with 


letting 


u  =  1  4-  to  cos  y  4-  Jii]  4-  J2u*  4  .  . 
y  =  0  -  wv0  4-  Jy  1  4-  J'y:  4-  . . . 

;  =  ;  0  4-  J 1 1  +  J  2 1 2  4-  .  . . 


An  algorithm  for  the  perturbation  procedure  1-: 


Let  n  =  1 


Svbstiivt  e  expression*  (ds  i~ /:]())  m/o  the  equations  of  motion  (2bi-(27) 


Equate  the  coefficients  of  Jr‘ 

Choose  the  arbitrary  constants  so  secular  terms  unll  not  arise. 
Solve  for  the  nth  order  solution 
Satisfy  the  inttial  conditions  (If)- (18) 

Iterate  on  a 


* 


The  calculations  were  carried  out  with  the  symbolic  manipulation  program  MACS1*  MA. 
In  this  paper  we  only  briefly  outline  these  calculations:  for  more  details  see  the  theses  of 
Sagova^  and  Snider. 

Beginning  by  substituting  equations  (28)  and  (30)  into  (26).  and  equating  the  term.- 
multiplied  by  J .  we  obtain 

~  —  —sc sin  2 ft - —  sin(t/  +  20)  -t - -  sin(y  —  20 )  (31  j 

do  2  '  2 

A  solution  to  this  equation  is 

?i  =  —  cos  2 0  +  — ^cos(y  +  20)  -f  — ~  cos(y  -  20)  +  A' ,  cos(2y  -  20)  4-  K2  (32) 
2  6  2 

The  last  two  terms  may  be  added  because  they  are  to  lowest  order  homogenous  solutions 
to  equation  (30).  The  term  multiplied  by  the  constant  A',  was  added  to  eliminate  secular 
terms  in  ? 2:  note  that  differentiating  this  term  with  respect  to  0  produces  terms  multiplied 
by  J.  from  equation  (20).  The  constant  Ah  was  added  to  satisfy  the  initial  condition  (17;. 
which  implies  that  ii(0u)  =  0  so 


A' 2  =  cos 20o  -  7—  cos(30o  -  -.’o)  -  — r"  cos(0o  +  ~0 )  -  h\  cos  2~0 

2  6  2 

Substituting  equations  (28 )- ( 30 )  and  (32)  into  (27).  and  equating  terms  multiplied  by  ./ 
yields 


(Pu\ 

~d(P 


3.' 


+  1/1  —  1  — 


+  (n  - 


5.s2 


—  +  1  j  +  -[(2  +  bt2)s2  -  2c2jcos2 0 


-  ».  2 

+  —  ( —  9.s2  +  8)  cos  2y  +  —  ( 1  Is2  —  6)  cos(  v  +  20)  +  — —  ( 3>2  —  2 )  cost 2 //  —  20 1  (33 1 

4  ‘  3  24 


0 


|(3s!  -  2)  -  Ml 

o  c 


cos(2 y  -  20) - — -  +  co  ^2^  +  4  -  5s2  j  cos  y  +  sin  y 


In  the  above  equation,  the  cos  y  and  siny  terms  would  produce  secular  terms  0siny  and 
Ocosy  in  uj.  The  choice  dy\/dO  =  5s7 /2  —  2  will  eliminate  these  possibilities.  Integrating 
yields 

y\  -  ^ —  2^  (0  -  0Q)  +  A’3[sin(2y  -  20)  +  sin  2ui0]  (34) 

The  term  multiplied  by  A'3  was  added  to  eliminate  secular  terms  in  u2-  The  constant  terms 
in  (34)  were  added  to  satisfy  the  initial  condition  y{0o)  =  0o  —  u?0. 

A  solution  to  lowest  order  of  equation  (33)  is  then 


ui  =  1  — y  +  T-  +  +  +  +  2e°l 


cos  20 


2  2 

+  ~(9*2  -  8)  cos2y  +  ^(-lls2  +  6)  cos(y  +  20)  +  ^(-3s2  +  2)  cos(2y  +  20)  (35) 

12  24  24 


o  2s h\ 


-2(3s2-2)- 
o  c 


cos(2 y  —  20)  —  -  I  -  +  I\4  cos {y  —  20) 


+/\5  cos  (y  —  Oq  4-  u,’o)  +  A  6  sin(y  —  Oq  +  u-'o) 


The  term  multiplied  by  K4  was  added  to  eliminate  secular  terms  in  u2.  The  terms  multiplied 
by  A'5  and  J\6  were  added  to  satisfy  the  initial  conditions  (1!)-(16). 

With  all  terms  in  place  to  deal  with  secular  terms,  the  calculations  are  continued  by 
substituting  equations  (28)-(30),  (32).  (34).  and  (35)  into  (26)  and  equating  terms  multiplied 
by  J 2: 

di2 

To 

We  have  for  brevity  only  indicated  on  the  right  side  of  equation  (36)  the  term  that  would 
produce  secular  terms  in  i2 .  Removal  of  this  term  by  making  its  coefficient  zero  determines 
A'j.  Equation  (36)  is  then  integrated  to  determine  ?2. 

Continuing  the  procedure  by  equating  the  terms  multiplied  by  J2  in  the  expansion  of 
equation  (27)  determines  y2.  A'3.  and  I\4.  Final  values  of  all  the  constants  are  listed  in 
Appendix  I. 


A'i  + 


scco(15s2  -  14) 
24  (5s2  —  4) 


sin(2y  -  20)  + 


(36) 
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Knowing  the  solution  for  i(6 j.  we  can  determine  Q(0)  by  integrating  equation  IS)  and 
applying  the  initial  condition  (18).  The  angle  6.  which  increases  continuously  from  an  initial 
value  #0<  may  be  related  to  the  time  t  by  num'okatly  integrating  (22). 

5  Solution 


Here  we  assemble  the  complete  solution: 


=  pol  1 1  +  e0  cos  y  +  J 


3s 


*>  °T  4/12 


1  —  — b  |  1  — 

to 


os2 


+  -y(-(2  -f  otl)s2  +  2cq)  cos  26 


+  r^(9s2  —  8)cos2y  +  -^(-lls2  +  6}cos(y  +  20)  +  r^(-3s2  +  2)cos(2y  + 
I Z  ZA  ZA 


2  6\ 


+  -^(3s2  -  2)  cos(2y  —  26) 


+ 


e0[lo(2  4-  e^s4  -14(4  4-  Cq)s2  4-  24]  sin  y  (os2  -  4)]  sin [6  4- 


k,’o] 


12(os2  -  4 ) 


+  ' 


CqS2(  15.«2  -  14)  sin  y  (os2  -  4);  sin  ;2~0  -  y  (os2  -  4)j  t^s2 


6(os2  -  4) 
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=  2,2 


+  ^“(3-s2  —  2)  cos (y  —  3^o  +  3~0 )  -  -—r  cos (y  -  5 0q  4-  3_0 ) 

e  3(  <;2 

+  y(3s2  -  2) cos(y  -  260  4-  2*c0) - cos (y  -  4 0O  4-  2-.’0) 

4  8 


cos(  y  -  6o  4-  3u.-0 ) 


(37) 


— -(s2  4-  1 )  cos (y  +  2-.-0)  4-  -((  —  2  4-  5c^)>2  —  2c2j  cos{y  4-  <?o  4-  ~'o! 

4  c 

+  -[(6  -f  oc.q)s2  —  4(1  +  Cq)]  cos(y  —  #o  +  >*-'o) 

+  —  [—(14  -f  5e2)s2  4-  2e„]  cos (y  —  3 60  4-  •‘•'o) 

2  2 
+  ^(9s2  —  4)  cos(y  4-  30o  —  «^'o)  4-  ~(  — "s2  4-  6)  cos(y  4-  9q  —  -2o) 

e2 

+  rr (~5s2  4-  4)  cos(y  -  60  -  -20) 
lb 

+  y(2s2  —  1 )  cos(y  +  20o)  +  ~(  — 3s2  +  1 )  cos(y  —  2#o)  4-  — (  —  3s2  4  2)  cos  y 
4  4  4 

+  e0.S2  COS($@  4  u20  )  4  yj-  cos(30o  -  «2o)  +  f-2  cos  20oj  j  4  p00(J2.  J36) 
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y  =  0  -  4jo  +  J  ( -  2  j  (0  -  0O) 


Jt l  f  (-75s6  +  260s4  -  296s2  +  112)  sin  y(5s2  — 4)  cos  !2- 

24(5s2  —  4)  {  (5s2  -  4)  "  ~ 

j~2 

+J0s2(-15s2  +  14)(15s2  -  13) cos 2-u’o 


JO 


( 5s2  -  4 ) ; 


M5? 


coS  ,  ,  .  2 


( 15s'  —  13)  cos(0o  -r  — 'c 


+  ^-(15s2  -  13)cos{30o-wo)  +  y(]5s2  -  13) cos 20o 
+  ^[5(9e2  +  34 )s4  +  4(9c2  -  34)s2  -  56e2]}  +  0(J2.  J30 ) 


(38) 


i  = 


to  +  scj|^  cos  2#  +  cos(i j  +  20) 

CqI  -  15s2  +  14)  sin  r4r  (5s2  -  4)]  sin  [2— -0  -  =r  (  5s2  —  4)1 

.  .  -- _ t* _ J _ [ _ + _ J 


1 2( 5s2  -  4) 


+  —  cos(y  —  20)  + 

cos  20o  ~  y  cos(3^o  ~  w0)  -  y  cos  (0O  +  wo)}  -f  0(J2.  J30) 


n  —  Qo  4-  cj 


0o  —  0  +  -  sin  20  -  e0  sin  y  w  -2  sin ( y  +  20)  — y  sin(y  —  20  i  -  ]■  sin  29 u 
2  6  2  2 


+eosin(0o  -  w0)  -  sin(30o  -  w0)  -  ^  sin( 0O  -4  w0) 
o  2 


f  2(  15s4  -  45s2  +  28 )  sin  y  (5s2  -  4 

)  cos  [2w0  -  y  (5s2  -4)1 

t  *  J 

l  (5s2 -4) 

+ - ^ - 

12(5s2  -  4)1 

+J0s2(los2  —  14)cos2w0}  -f  cJ20 


e0s 


-e0s  cos(0o  +  w0) - —  cos(30o  -  w0) 


3 


-s2  cos  200  +  ~~(7s2  -  4)  -f  —  (-S2  +  6) 

1 1 


k  +  O(J2.J30) 


(40) 


1  r9 

t  =  (o  +  7"  /  c2  (  1  +  J 
0  o 


(-3s2 +  2) 


cos  20  +  to(s2  -  1 ) 
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COS (y  -  2 8) 


e0(~4.s2  +  3)  e0(-2s2  +  l) 

cos  y  + - - - COS  (y  +  20)  + - - - 


eQ.<i2(15s2  —  14) sin  y(5s2  —  4)  sin  2-v0 

_____ 


f  (5s2 -4)! 


+5  —  1  +  —  cos  2@o  H"  ~z —  cos (3^o  —  ^’o)  *4 — ~ —  cos($o  4-  ^o)  4*  —0[  J  ,  J  0 ) 

1  6  l  J  J  fto 

In  obtaining  the  equations  (37)— (41 ).  use  has  been  made  of  trigonometric  formulas 
to  simplify  terms  containing  the  factor  5s2  —  4  in  the  denominator.  In  the  form  given, 
these  terms  can  clearly  be  seen  to  approach  a  finite  limit  at  the  ‘'critical  inclination" 
i0  —  sin-1  ^4/5  =  63°26'  or  116°34'.  Hence  the  solution  is  actually  valid  for  all  values 
of  i'q.  If  | ?o  —  sin-1  \J^/o\  <  J.  the  formulas  ( 37 )- ( 4 1 )  can  still  be  used  by  letting  os2  —  4  —  J. 
or  the  limiting  forms  for  ?0  — +  sin-1  ^4/5  can  be  used. 

To  check  the  solution,  we  can  see  if  the  specific  mechanical  energy  (18)  of  the  satellite 
remains  constant.  Substitution  of  the  solution  (36)-(37)  into  equation  (10)  plus  (11)  yields 


r  +  r  =  - 


GM ( 1  -  t2 )  GUJ2R2(  1  -  3 sin2  30 )  GM  "  ;2 


+ - 0{JZ) 

Po 


The  right  side  is  easily  recognized  as  the  value  of  the  specific  mechanical  energy  at  the  initial 
time  t0- 

As  a  further  check  on  the  solution,  we  can  see  if  it  reduces  to  our  previous  results  for 
equatorial  and  polar  orbits,  obtained  by  completely  separate  derivations  (Danielson  and 
Snider,  1989).  Setting  to  =  0  and  using  the  independent  variable  a  measured  from  the  line 
07,  we  find  that  equations  (37)— (41 )  reduce  to  equations  (18 )— (22 )  of  our  previous  paper. 
Setting  t0  =  7t/2  and  using  the  expansion  cos (y  +  Jk)  cos y  -  Jks'my.  we  find  that 
equations  (37)— (41)  reduce  to  equations  (38 )— (4 1 )  of  our  previous  paper. 

Comparing  the  terms  in  the  0-symbols,  we  see  that  the  relative  error  in  equation  (41) 
may  be  greater  than  that  of  equations  (37)-(40).  Since  the  underlined  terms  in  equations 
(37)-(40)  are  of  this  same  order  of  magnitude,  we  can  drop  the  underlined  terms  except 
when  (37)— (38)  are  used  to  calculate  r  in  equation  (41).  The  relative  error  of  our  solution 
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will  then  still  be  of  order  (0  —  Oq)J2. 

If  we  retain  only  the  two-body  solution,  the  relative  error  terms  will  be  of  the  ordei 
(6  —  60)J.  Here  the  error  in  our  solution,  as  compared  to  the  exact  solution  of  the  equation- 
of  motion,  should  be  of  the  order  J  times  the  error  in  the  two-body  solution  (for  an  Earth 
satellite  J  <  .0015). 

6  Comparison  of  Perturbation,  Two-Body,  Numerical,  and  Mea¬ 
sured  Solutions 


In  this  section  we  compare  the  preceding  perturbation  solution,  the  two-body  solution,  a 
completely  numerical  solution  of  the  differential  equations,  and  actual  measured  satellite 
data;  for  more  comparisons  see  the  thesis  of  Krambeck.  The  difference  between  the  position 
vector  r  determined  by  the  numerical  integration  code  or  measured  data  and  the  position 
vector  rref  calculated  from  our  perturbation  solution  or  the  two-body  solution  is  the  error 
Ar. 


Ar  =  r  -  rr„f 


If  the  errors  (Sr.  SO.  At.  AD)  in  the  orbital  parameters  (r.  0.  i.  ft)  are  small,  we  can  estimate 
Ar  from  equation  (4)  and  the  linear  approximation 


.  dr  dr  dr  dr 

Ar  %  — Ar  +  —A 6  +  —At  +  —AH 

dr  d0  di  dP 


ir>) 


It  is  customary  to  decompose  Ar  into  components  (6i<62.  S3)  along  the  moving  triad  (B, .  B2.  B . 


Ar  —  ^jB;  -+■  62B2  -f-  63B3 


The  component  61  is  called  the  radial  error.  82  is  the  down  track  error,  and  f>3  is  the  cros- 
track  error.  Applying  (42)  to  equation  (4).  and  expressing  the  base  vectors  (b].b2.b3)  in 
terms  of  (Bj.B2.B3).  we  obtain  the  following  approximations: 

61  Ar  .  82  %  r(  A#  -1-  cos  ? SP) .  m  r(sin  OSi  —  cos  0  sin  ?Af>)  (43i 
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We  obtained  the  numerical  integration  code  UTOPIA  from  the  ('dorado  Center  for 
Astrodynamics  Research  located  on  the  campus  of  the  University  of  Colorado.  The  code 
was  specialized  to  the  differential  equations  used  in  this  paper.  We  compared  the  solutions 
for  an  earth  satellite  with  the  following  initial  conditions: 

r0  =  7,386.18  km 
t0  =  .003991 
60  =  104.05° 

«;o  =  224.38° 
i'o  =  90.03° 
n0  =  322.63° 

to  =  0 

These  initial  conditions  represent  an  essentially  polar  orbit  at  an  altitude  of  approximately 
1000  kilometers  and  period  about  l|  hours.  For  this  satellite  the  perturbation  and  numerical 
orbits  match  extremely  well  while  the  two-body  orbit  is  grossly  erroneous.  The  magnitude  of 
the  error  in  r  is  shown  in  Figure  2.  Note  that  the  relative  error  in  our  perturbation  solution 
is  2.8 J2{6  —  0o).  and  that  this  error  is  1.1J  times  the  error  in  the  two-body  solution. 

We  obtained  measured  satellite  data  from  the  First  Satellite  Control  Squadron  located 
at  Falcon  Air  Force  Base.  Colorado.  A  near  earth  satellite  processed  the  following  initial 
conditions: 


r0  = 

7.  776.58  km 

to  = 

.0003071 

II 

O 

149.14° 

U,'o  — 

9.57° 

io  = 

98.81° 
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n0  =  37.10° 


to  =  0000Z  26  July  1990 

Again,  the  perturbation  orbit  is  far  superior  to  the  two-body  orbit.  The  radial,  down  track, 
and  cross  track  errors  (6ly  63,63)  are  shown  in  Figure  3.  Note  that  although  the  perturbation 
solution  produces  only  a  small  iropro%'ement  in  the  radial  error,  this  error  is  negligible  in 
comparison  to  the  down  track  error. 

7  Conclusions 

Our  solution  embodies  the  principles  outlined  in  the  introduction.  The  relative  error  of  our 
solution  is  of  order  (6  —  0o)J2 ,  which  is  a  factor  of  J  times  the  relative  error  of  the  two-body 
solution;  our  solution  loses  its  validity  after  an  angular  change  (6  -  60)  of  order  l/./2.  which 
is  a  factor  of  j  longer  than  the  interval  of  validity  of  the  two-body  solution.  Secondly,  our 
solution  is  in  terms  of  classical  orbital  elements:  no  transformation  to  an  alternative  non¬ 
physical  set  of  elements  is  required.  Finally,  our  solution  is  free  of  singularities  for  all  value- 
of  the  initial  orbital  parameters,  including  elliptic,  parabolic,  and  hyperbolic  orbits. 

Our  formulas  should  agree  closely  with  satellite  orbits  whose  dominant  perturbation  i- 
the  planet's  oblateness.  Of  course,  the  effects  of  higher-order  terms  in  these  expansion-, 
higher-order  terms  in  the  planet's  potential,  and  of  other  perturbation  forces  may  also  be 
important.  The  formulas  will  have  to  be  amended  to  include  these  additional  effects. 

APPENDIX  I:  Values  of  the  Constants  K^-Kq 


A’,  = 


csCq(  —  15s2  +14] 
24 (os2  -  4) 


AN  = 


sc£o  ,nn  ,  scr-o  ,rt  ,  ,  ,  c.se2(15.s2  -  14) 

-  —  COS  200 - —  COS(30q  -  a?0) - —  COS ( +  w'o)  +  — -  r - - -  COS  2w'0 

2  6  2  24(o52  —  4) 


AN  = 


cg(-75^  +  260s4  -  296s2  +112) 
48 (os2  -  4 )2 
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,,  [15(e2  +  2)54-14(e2-M)52  +  24; 

^4  ^0  O  <  i  *  - 2  1  \ 

24 —  4) 

2  2 
=  y§(-9*2  +  8)  cos(20o  -  2u-'0)  +  ^(3s2  -  2) cos(4#o  -  2-;0) 

e  t2 

—  (eos7  +  A4)  cos(0o  4-  ^'0)  ■+•  — (s2  ~  2)  cos(30o  —  >*0}  4-  —(—3 6*  +  2)  cos  2*^o 

b  b 

—  —  (5(2  —  Cq)s2  4-  2cq]  cos 20o  +  pU'^eo  4-  18)s2  —  (cq  4-  1 ) 

2  2 

4-~r(6.s2  —  5)  sin(2#0  —  2-Jo)  4-  —( —  3.s2  4-  1 )  sin(4#0  —  2—’o) 
o  12 

4'^(eo(—  4-1)4-  2A4]  sin(0o  4-  -•‘0)  4-  —  (352  —  2)  sin($o  —  ~o) 

e  e2  1 

4-  —  (  —  75 2  +  2)  sin (3#0  —  >*'o)  4-  -— ( —  -s2  4-  1 )  sin  2u,'0  4-  —  [— (ocq  -f  2 ) 5 2  4-  2e2]  sin  2$o 
b  4  o 


APPENDIX  II:  Rigorous  Bounds  on  the  Orbit 

It  follows  from  ( 1 0 )-( 1 2)  that 

_  1  fdrV  r  d7r  GM  GUJ2R2iy  o.2. 

T+l  =2(*)  +2^-^-  +  ^C-|1-3s'" 

This  can  be  rewritten  in  the  form 

~  r7(d-^\  =  4(7  -f  \~)r  -f  26’ -U  4-  —  7i2/?-  (3sin2i-  1) 

dr  \  dt )  r* 

from  whence  it  follows  that 

±1>(±Y)  <4(r+l>  +  2C.W  +  gjgg 

dr  \dt )  rz 

Integrating  from  r(t0)  to  r(t)  yields 

r*  MV  <  2,r  4  w  +  2CA/r  -  -  *;  +  HIM. 

\  dt )  r  Poro 


It  follows  that 


0  <  2(T  +  \  ')r2  +  2GMr  -  h20{\ 


3J2R2. 
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When  T  +  V  <  0,  the  quadratic  polynomial  on  the  right  side  of  (44)  has  the  roots  (exact 
values  can  be  found  from  the  quadratic  formula) 


rmin  —  - - [1  +  0(^2)]  •  Cmax  =  - [1  +  0(J2l\ 

1  T  to  l  —  to 

Hence  a  satellite  having  negative  mechanical  energy  remains  for  all  time  within  the  spherical 
annulus  r^n  <  r  <  r^.  Since  the  position  vector  is  bounded,  we  can  invoke  the  recurrence 
theorem;  i.e.,  the  satellite  will  come  as  close  as  desired  to  its  initial  position  in  a  sufficient!} 
long  period  of  time  (as  shown  by  Poincare).  Furthermore,  we  are  guaranteed  of  the  validity 
of  supressing  secular  terms  to  describe  the  orbit  via  perturbation  analysis. 
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